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The theory of nonhnear diffraction of intensive light beams propagating through photorefractive 
media is developed. Diffraction occurs on a reflecting wire embedded in the nonlinear medium at 
relatively small angle with respect to the direction of the beam propagation. It is shown that this 
process is analogous to the generation of waves by a flow of a superfluid past an obstacle. The 
"equation of state" of such a superfluid is determined by the nonlinear properties of the medium. 
On the basis of this hydrodynamic analogy, the notion of the "Mach number" is introduced where 
the transverse component of the wave vector plays the role of the fluid velocity. It is found that 
the Mach cone separates two regions of the diffraction pattern: inside the Mach cone oblique dark 
solitons are generated and outside the Mach cone the region of "ship waves" is situated. Analytical 
theory of "ship waves" is developed and two-dimensional dark soliton solutions of the equation 
describing the beam propagation are found. Stability of dark solitons with respect to their decay 
into vortices is studied and it is shown that they are stable for large enough values of the Mach 
number. 

PACS numbers: 42.65.-k, 42.65.Hw, 42.65.Tg 



I. INTRODUCTION 



An analogy between propagation of light beams in nonlinear media and superfluid flow is well known and quite sug- 
gestive. Formally, it is based on a mathematical similarity of the equations for electromagnetic field evolution of light 
beams in paraxial approximation and Gross-Pitaevskii equations for superfluid motion of Bose-Einstein condensates 
of dilute gases. Accordingly, such nonlinear structures as bright or dark solitons and vortices have been thoroughly 
studied both in optics and superfluid dynamics (see, e.g., [Hll]). These structures arise as a results of interplay of non- 
linear and dispersive properties of the medium under consideration. One more example of such a structure is provided 
by so-called dispersive shocks which replace a notion of usual dissipative shocks in compressive fluid dynamics in case 
when dissipation can be neglected compared with dispersive effects. As a result, a thin layer with strong dissipation 
within unfolds into a region with fast oscillations, which can be represented as a modulated nonlinear periodic wave 
(a "soliton lattice"). The notion of dispersive shocks arose first in water wave physics (the theory of undular bores 
on rivers) [3j and plasma physics (coUisionless shock waves) |5] , then generality of this phenomenon was realized and 
(based on the Whitham theory of modulations of nonlinear waves) mathematical methods for their description 
were developed [6]-[T2]. 

Realization of Bose-Einstein condensate of dilute cold gases [13j [HI [15] and study of its dynamics has naturally 
led to the theoretical and experimental studies of dispersive shocks in this new medium [IB]-[5D]. Corresponding 
optical counterpart of dispersive shocks suggested by the above mentioned analogy between beam optics and superfluid 
dynamics was realized experimentally in [211 [22 [231 Ell and the theory of such optical dispersive shocks was developed 
in 125). 
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FIG. 1: A sketch of formation of nonlinear diffraction pattern in propagation of a light beam through photorefractive medium 
with embedded reflecting wire. 



In dissipative fluid dynamics with negligible dispersion shocks can also be generated by a supersonic flow of the 
fluid past an obstacle. Such shocks have the form of a sharp stationary jump of the fluid parameters across certain 
lines inclined with respect to the flow direction. For shocks of small intensity these lines lie along the so-called "Mach 
cones" (see, e.g. |26]). In dispersive fluid dynamics these oblique shocks unfold into "fans" of spatial solitons spreading 
downstream from the obstacle ^27j. The theory of such oblique dispersive shocks was developed in [28J for the case 
of weakly dispersive media when the flow past a slender body is asymptotically described by the Korteweg-de Vries 
equation along the Mach lines. 

Dynamics of a Bose-Einstein condensate is described by the Gross-Pit aevskii equation and the theory was extended 
to this case in pOj. If the obstacle is small enough, then the shock consists of a single oblique dark soliton. The theory 
of oblique dark solitons was developed in [301 131] • It is important to note that such oblique solitons are located inside 
the Mach cone with the Mach number deflned as a ratio of the flow velocity to the sound speed calculated at inflnite 
wavelength. The so-called "ship waves" arising as stationary dispersive wave packets of Bogoliubov excitations are 
located outside the Mach cone. Apparently, they were observed in the experiment |31] and their theory was developed 
in [331 1311 ■ The analogy between beam optics and superfluid dynamics suggests that similar effects would exist in the 
optical context where they take the form of diffraction wave patterns in light beams propagating through a nonlinear 
medium. Although such structures were observed in some experiments (see, e.g., |35p. they have not been studied 
systematically yet. In this paper, we shall consider a typical simple situation of nonlinear diffraction of light which 
can be considered as an optical counterpart of generation of spatial dispersive shocks and "ship waves" in the flow 
of Bose-Einstein condensate past an obstacle. To be definite, we consider a light beam propagating through a bulk 
self-defocusing nonlinear refractive medium with a thin wire (a "needle") inserted in it; see Fig. 1. Direction of the 
light beam is tilted with respect to the wire that is there exists a "flow" of light "past an obstacle". As a result, at 
the output plane of the medium a diffraction pattern is formed consisting of oblique dark solitons and "ship waves" . 
We shall give here analytical and numerical treatment of this phenomenon and obtain main characteristics of the 
diffraction pattern. 



II. MAIN EQUATIONS AND GENERAL FORM OF THE DIFFRACTION PATTERN 



Propagation of stationary beams is described by the equation 

■ dip 1 « , fco 

OZ ZKq Uq 
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where ip is envelope field strength of electromagnetic wave with wave number fco = 27rno/A, z is the coordinate along 
the beam, x,y are transverse coordinates, r = (x, y), A_l = d^/d^x + d^/d^y is transverse Laplacian, no is a linear 
refractive index, V{v) represents a "potential" of an obstacle (e.g. a reflecting wire) at which diffraction occurs, and 
in a photo-refractive medium we have 

5n=-\nlr:,^Ep^^, (2) 

^ P + Pd 

where Ep is applied electric field, r^^ electro-optical index, p — and pd is the saturation parameter. 
For mathematical convenience, we introduce non-dimensional variables 

1 ;„ 2 t:. f Pc \ ~ / 1 „ t:. ( Pc\ ^ ~ / 1 „ t:. ( Pc 



-kuQrssEp y—j z, X kno^ -rssEp y—jx, y = kno]^ -rj^Ep ^—jy, ip ^ ^t};, (3) 

where pc is a characteristic value of optical intensity (its concrete definition depends on the problem under consid- 
eration; for instance, it can be the background intensity), so that Eq. ([T]) takes the form of generalized nonlinear 
Schrodinger (GNLS) equation 

.g + ^A.^- J^^ + nr)V. = 0, (4) 

where 7 = pd Pd, ^(r) is represented in non-dimensional units, and tildes are omitted for convenience of the notation. 
In fact, our approach can be applied to other forms of the nonlinear term provided it corresponds to self-defocusing 
light beams. Therefore we shall also use the general form of the equation 

z^-f iA^^-/(|^|2)^ + l-(r)^ = 0, (5) 
oz 2 

where /(IV'P) > Oj I^i particular, for photorefractive medium, 

/(p)=p/(l-f7P)- (6) 

If saturation effect is negligibly small (7|'!/'P ^ l), then Eq. (|4| reduces to the standard cubic nonlinear Schrodinger 
(NLS) equation 

i'^ + \^^^-\iP\^^ + V{v)^ = {). (7) 

If the phase of -0 is a single-valued function, then it is convenient to represent the above NLS equations in a fluid 
dynamics type form by means of the substitution 



V'(r, z) — ^/p exp J u(r, z) ■ dr 
so that they are transformed into 

Pz+W±{pu) ^0 

u, + (uVj.)u + V^/(p)-Vnr)-Vi ^' 



(8) 



A^P _ (V^l ^ (9) 
4p 8p2 J 

In the hydrodynamic interpretation the light intensity p has a meaning of a density of a "fluid" and Eq. ([6| can be 
viewed as an "equation of state" for such a fluid. The function u(r, z) is a local value of the wave vector component 
transverse to the direction of the light beam; in hydrodynamic representation it has a meaning of the "flow velocity" . 
The variable z plays the role of time so it is natural to describe the deformations of the light beam in evolutionary 
terms. We note that substitution |8| rules out vorticity so that system ^ actually represent a restriction of the 
multi-dimensional GNLS equation (Isl to potential "flows" . 

We shall consider propagation of a tilted light beam with uniform input intensity, that is at z = it has the initial 
form 



^{r, 0) = exp{iUx), 



(10) 
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FIG. 2: Evolution of the diffraction pattern at tlie output plane as a function of the length z of the photorefractive medium. 
The patterns are obtained by numerical solution of Eq. Q with ^(r) corresponding to an ideally reflecting wire with unit 
radius for 7 = 0.2, t/ = 2, and (a) z = 20, (b) 2 = 40, (c) z = 60. 



that is we suppose that the background intensity is equal to unity; IJ represents the x-component of the wave vector 
due to tilting of the light beam. The problem is to describe the wave pattern at the output value of z. 

To clarify a general picture of the diffraction pattern, we have solved Eq. Q numerically for the initial wave 
function -0 given by Eq. (10 1 with U = 2 and the boundary condition of vanishing ■0 at the surface r = 1 of the 



obstacle located at a; = 0, y = 0. As we see, the diffraction pattern consists of two different parts separated by the 
Mach (or Cherenkov) cone which is defined as lines drawn at angle 6 with respect to the direction of the flow (x axis) 
with 



sin 9 



1 



M 



U 



where the sound velocity corresponds to the dispersionless limit of Eqs. ^ that is {yp/ p = \7f(p)) 



Po 



Po 



which in the photorefractive case with pq = I yields 

1 



1 + 7 



and M=f/(l + 7). 



(11) 



(12) 



(13) 



Outside the Mach cone, there is a stationary wave pattern created by interference of linear (far enough from the 
obstacle) waves. Inside the Mach cone there are two oblique dark solitons situated symmetrically with respect to the 
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FIG. 3: Distribution of the phase in the difTraction pattern at the output plane of the photorefractive medium. The pattern 
corresponds to 7 = 0.2, [7 = 2, and z = 60. 



direction of the "flow" . These oblique solitons decay at the end points into vortices but closer to the obstacle they 
are described by a potential flow with jump of phase across them as it is demonstrated in Fig. 3. 

Our task now is to develop analytical theory for these two regions of the diffraction pattern and to compare it with 
numerical simulations. We shall start with the "ship waves" pattern located outside the Mach cone. 



III. DIFFRACTION PATTERN OUTSIDE THE MACH CONE 



If the size of the obstacle is much less than the wavelength of the pattern, then we can consider it as a point-like 
one and take the obstacle potential in the form 

V{x) = V^8(v). (14) 

Far enough from the obstacle, the amplitude of the wave pattern is small compared with the background intensity of 
the light beam. Hence, the wave pattern can be calculated by means of perturbation theory |36j . 

If we neglect the influence of the obstacle, then the function of a uniform light beam with the intensity po 
depends on z in the reference frame with U = as "0 oc exp(— i/(po)-2)- We exclude this dependence by introducing 
the substitution = ^ ■ exp{—if{po)z) so that ^' satisfles the equation 

*vI/, + iAV'+[/(po)-/(|*n]*-0. (15) 

In the same reference frame the obstacle moves with the velocity — U and generates diffraction waves which in the 
linear approximation are described by a small correction S'i' to the unperturbed wave function: \E' « y/pQ + S"^. Hence 
S'i' satisfies the equation 

iS^i, + iAiJ^- - €^(5^- + (5**) - Voy/p^5{r + Uz) = (16) 

where we have added the potential of the obstacle due to which linear waves are generated. In the stationary case, 
which we are interested in, the wave pattern moves with the obstacle, that is in the reference frame attached to the 
reflecting wire we have = \E'(r + XJz) and 

^(5^'(r + Uz) = (UV)(5^'(r + Uz). 
dz 

Introducing r' = r + Uz and omitting primes, we arrive at the equation 

i(UV)(5* + iAJ^- - c2(5«' + (5^-*) - Vo^poSir) = , (17) 
describing stationary diffraction pattern generated by the beam. 



FIG. 4: Coordinates defining a radius vector r and a wave vector k, normal to the wave front shown schematically by a curved 
line. 



Equation (17 1 can be solved by the Fourier method. We introduce the Fourier transform of the wave function 



(2^)2 



and obtain 

- (kU + e/2 + c2)<5*k - c^J^lk = ■K)Va^. 
Another equation is obtained by means of substitution k —t — k and complex conjugation: 

- cl6^\. + (kU - e/2 - cl)5^*_]^ = FoVpo- 



Solution of Eqs. ( |19|20| reads 



fcV2 - kU 
(kU)2-fc2(c2 + fc2/4)' 



(18) 

(19) 
(20) 

(21) 



Since 



we arrive at the following expression for the intensity perturbation in the output diffraction wave pattern created by 
propagation of light past a reflecting wire: 



5p = Vopo 



2„ikr 



cPk 



(kU)2 - fc2(c2 + fc2/4) + iQ (27r)2 



(22) 



where we have introduced an infinitesimal positive imaginary term +iO corresponding to the radiation condition for 
outgoing waves. 

Now we introduce polar coordinates (see Fig. 4) defining the components of the vectors r and k as 



X — r cos X, y = r sm x; 
kx = — fccosTy, ky — ksinr]. 



Simple transformation casts Eq. ( 22 ) to the form 



. VoPo 
Op^ 



TT JO 



(23) 



(24) 
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where 



fco = 2cs \/ Iv'P cos^ r/ — 1 = Csk{rf). 



(25) 



We can represent the integral ( 24 ) as a sum 



37r/2 



-Tr/2 



7r/2 



(i?7 



-7r/2 



37r/2 



dry 



7r/2 



and, noticing that the second term after substitution if ~ rj — n becomes equal to a complex conjugate of the first 
one, we rewrite it as 



6p = ^Re 



-7r/2 Jq 



To perform integration over k, we notice that the integrand function has a pole in the first quadrant, 



k = y^feg + iO = ko + iO, 



(26) 



(27) 



which gives the main contribution into the integral for cos(x + v) < 0. Indeed, taking a closed contour along the 
positive real axis of k with added quarter of the circle, which gives no contribution into the integral, and a path along 
positive imaginary axis which contribution 



k'^ + fcg 



(28) 



is decreasing with r much faster than the contribution of the pole (which is proportional to r see below), we 
obtain 



dp = im 

TT 



tt/2 



n/2 



-ikr cos(x+7/) 



drj, 



(29) 



where k is determined by the equation (251 (index "0" is omitted here). 
If the phase kr np, where 

ip{rj) = k{i]) cos{x + v), 



(30) 



is large enough, the integral (29 1 can be evaluated by the standard method of stationary phase. This condition is 
fulfilled far enough from the obstacle r — > cxd provided \k(ri) cos(x + t?)] 3> 1/r. The equation which determines the 
point of the stationary phase dip/drj = gives relationships for the angles (see Fig. 4) 



2U^ 2AP (c2 + A;V2)tanr7 (1 + fcV2) tan 77 
tan = -— sm2r] = — — sm2ry, tanx=77T5 , „ = . 



(31) 



Taking into account equation ( 25 ) , we find 



cos p 



2[(M2-2)fc2+4(M2-l)]i/2 



With account of (31 1, we get the expression for the second derivative of the phase 



d^Lp „cosp 



fc3 



[(Af2 - 2)k^ + 6{M^ - 1) 



As a result, the expression for the condensate density ( |29| takes the form 

Sp = VqPq' 



' 2k [(M2-2)fc2 + 4(M2-l)] i/4 
7rr [(M2 - 2)fc2 + 6(M2 - 1)] 



1)11/2 V 



Cgkr COS p 



(32) 



(33) 



(34) 
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where 



k 2^1/2 cos2?7- 1. 



As we see from Eq. ( 34 1 , the Hnear waves exist only in the region 

— arccos(l/M) < ij < arccos(l/M) 

outside the Mach cone. 



(35) 



(36) 



With the help of Eqs. (31 ) one can find the shape of the hnes of constant phase (e.g. wave crests) <& — fcr cos fi in 
a parametric form 



X = r cos X ■ 
y = r sin x = 



4$ 
4$ 



c 



COS 77(1 - cos 277), 
sin r7(2M^ cos^ 77 - 1). 



(37) 



Small values of rj correspond to waves in front of the obstacle. In this case we have 



y '■ 



(2AP - i)$ 
'2c,VM2-1 ^4c,(M2- 1)3/2^ 

(2M2 - 1)$ 



(38) 



2c,(Af2- 1)3/2'^' 
that is the lines of stationary phase take parabolic form 



C,,(Af2 - 1)3/2^2 



2c,VM2 - 1 (2Af2 - 1)$ 
The limiting values i] — zt arccos (l/Af) correspond to the lines 



X 

y 



(39) 



(40) 



i.e. far from the obstacle the lines approach to the straight lines parallel to those forming the Mach cone (111. 
Predictions of the analytical theory are compared with the numerically calculated wave pattern in Fig. 5 and excellent 
agreement is found. 

In the region in front of the obstacle where ?/ = 0, x < 0, the perturbations of the light intensity take the simplest 
form. Here we have 



k = 2csVm2 - 1, 



i.e. the wave length A = 27r/fc is constant and 

' (M2- 1)1/2 



5p 2Vf)PQ\ 



7r(2Af2 + l)|x| 



IS {-2cs\f- 



M2 - 1 ; 



^) , y = o, x<Q. 



(41) 



(42) 



The plot illustrating this dependence is shown in Fig. 6. As we see, approximate formula Eq. (34 1 is accurate enough 
almost everywhere except the small vicinity of the obstacle. 

As was indicated above, the method of stationary phase used for the derivation of (34 (requires the condition 
|/c(?7) cos(x + 77) I ^ 1. According to (251 we have /c ^ at the Mach cone and the necessary condition is not fulfilled. 
To find a wave pattern near the Mach cone one should return to the investigation of the integral ( 22 1 and introduce 
new coordinates along the Mach cone (^) and normal to it (r) (i.e., they are rotated to the angle 9 around the origin): 

X = ^cos 6' — r sin6', y = ^sin6' + rcos^. (43) 



In new coordinates equation (22 1 takes the form 



5p = Vopo 



dkfdkr 



{k^U COS0- krUsiney - k^{c^, + k^/A) + iO (27r)2 



(44) 
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FIG. 5: Numerically calculated wave pattern corresponding to diffraction of a light beam on the obstacle embedded into a 
photorefractive medium. The plot corresponds to 7 = 0.2, U = 2, and the radius of the reflecting wire to r = 1. Dashed line 
corresponds to linear analytical theory, Eq. (371, for the line of constant phase; it is shifted to the left to two units of length 



from the center of the obstacle due to its finite size in numerical simulations and better fitting to numerics. 




FIG. 6: Profile of intensity in front of the the obstacle for x < 0, y — and choice of the parameters 7 — 0.2, U = 2, Vo — 2.6. 
Solid line corresponds to Eq. (|42|| and dashed line to numerical solution of Eqs. (|5|l 



Far from the obstacle, near the Mach cone, the dependence of the wave pattern on the ^-coordinate is much slower 
than dependence on the r-coordinate; besides that one has <C 1 here. Main contribution into the integral over /cj 
is due to the pole which position is determined by the equations 



(k^U cose- krUsinef - fc2(l + k^/A) = 0, 



Their approximate solution for fc^ ^ /c,- <C 1 is given by 

h = - 



~\~ — - k . 



(45) 



(46) 
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where we have taken into account Eq. (13 1. Integration over yields 



Sp^ 



Vo/Oo d 



- 1 dr 

and with account of the integral representation of the Airy function 

poo 



COS I , — KrT dkr 



M{z) 



cos + zhi^ dn 



we obtain the following expression for the density oscillations in the vicinity of the Mach cone: 



2VoPo 



(M2 - 1)1/6(3^)2/3 



Ai' 



2(Af2- 1)1/6 



(47) 



(48) 



(49) 



where Ai' denotes the first derivative of the Airy function with respect to its argument. Returning to x and y 
coordinates, we get 



dp- 



'^VoPo 

(M2 - l)i/6[3(2:cos6i + ysin6i)]2/3 



Ai' 



2(M2- 1)1/6 

— : — TTTT-^ [-X sin fc* + y cos ( 

[3(a::cos6' + 2;sm6')]i/3 ^ 



(50) 



where sin6l = 1/M, cos 9= ^/AP - 1 /M. 

The above formulae allow one to derive expressions for the dependence of intensity on y coordinate for fixed value 
of X which may be convenient for comparison with the experiment and numerical simulations. Far enough from the 
Mach cone when Eq. ( 34 ) can be applied we find dependence of x on y from the equation 



y (l + fc2/2)tann 
— — tan ri — — . 



then 



r(r/) 



sinx(?7) ' 



(51) 



(52) 



where k{rf) and p.{ri) are defined by Eqs. (35) and (32l. In the limit y x we have x ^ '''/2, hence denominator in 
the rhs of Eq. (51 1 vanishes and 



k ^ ^/2{NP - 1) for y > X. 



Comparison with Eq. (35 1 gives the limiting value of rj, 



COST] ' 



A/M2 + 1 

V2M 



Substitution of these values of the parameters into Eq. (|34|) yields 



2M 



n{AP - l)y 



cos [(M - 1 /M) y - 7r/4] for y > x. 



(53) 



(54) 



(55) 



The profile of the wave in the vicinity of the Mach cone is shown in Fig. 7. As we see, Eq. (34 1 reproduces the density 
profile very well almost everywhere except for a closest vicinity of the Mach cone and inside it where the density 
perturbation decays exponentially according to the behavior of the Airy function in Eq. ( 50 1 . 



IV. OBLIQUE DARK SOLITON 



Far enough from the obstacle where vorticity is equal to zero and the light "flow" can be considered as potential, 
we can use the hydrodynamic representation of equations of light beam evolution. Here the potential of the obstacle 
can be neglected (in case of a reflecting wire it obviously vanishes beyond the surface of the wire, that is the obstacle 
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FIG. 7: Wave pattern near the Mach cone. Solid line corresponds to Eq. (1341) and dashed line to Eq. ||50| 



is represented by an infinite cylindrical barrier) and for large enough z the soliton is close to its stationary state. The 
profiles of intensity p and "velocities" w, v can be found analytically as a solution of stationary equations 



+ {.pv)y = 0, 



(56) 



and 



UUx + VUy + 



UVx + VVy + 



1 + 7P/^ 

p 

1 + 7P 



Px ~^ Py Pxx Pyy 



pI + pI 



4p 



= 0, 
= 0, 



(57) 



with boundary conditions (in this Section we assume po = 1) 

p=l, u = U, v = at |x| 



(58) 



To simplify calculations, it is convenient to notice that one of equations ( 57 1 can be replaced by the condition of zero 
vorticity 



Uy - = 

which is fulfilled for the potential fiow in the soliton solution. 
We look for the solution in the form 



p — p{9), u — u(0), v — v{9), where 9 — x — ay. 



(59) 



(60) 



The parameter a determines a slope of the oblique soliton in the x, y plane. Then equations ( 56 1 and ( 59 ) with account 



of conditions (|58| give after simple calculation the expressions for the components of the "fiow velocity" in terms of 

C/(l + a2p) aU{\-p) 



the light intensity 



{l + a?)p 



(l + a2)p 



Substitution of these expressions into any equation ( 57 1 and integration of the resulting equation yields 

1 + 7P V 2 1 + 7 



1 



(l + a2)2(p'2_2pp") + (l + a') ^' 



= 



(61) 



(62) 
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where an integration constant is chosen in accordance with the conditions ( 58 1 . This equation can be integrated once 
more to give 

2 



(l + a2)p 



7 



ln(l+7/9)- 



1 



(1 + 7)7 



P 



7 



■ln(l + 7) - 



1 



7(1 + 7) 



2 



(63) 



where (58 1 is also taken into account. For a given "Mach number" M = (l + 7)[/ the soHton solution depends on the 
slope parameter a alone. 

Now we notice that expressions for the flow velocity field (61) in terms of intensity p do not depend on the 
nonlinear properties of the medium but are determined completely by the "continuity" equation and the condition 
(59 1 of potentiality of the flow. Therefore we can change the "reference frame" in such a way that the transversal 
velocity (wave vector u) is equal to zero at \9\ — > 00. This means that we rotate the reference frame to the angle 
(j) = arctana and pass to the frame "moving" with "velocity" (U cos(f>, t/sin^) as z increases, which means the change 
of coordinates 



X - 
V 



XCOS ( 

- X sin ( 



ysm(; 
j/cos< 



Ucos ( 
- U sin ( 



Correspondingly, the "velocity" field transforms as 

u = {u — U) cos 



— V sm <p, 
{u — U) sin (f>+ V cos 0. 



Substitution of (61 1 gives 



where we have introduced the parameter 



- 1 



V = 0. 



U 



(64) 



(65) 



(66) 



(67) 



In new variables the velocity field does not have a component along y coordinate. The variable 9 takes the form 
9 = \/l + a'^{x + cz) and hence the intensity p does not depend on y coordinate. Thus, in new coordinate system we 
have a ID dark soliton moving with velocity c in negative direction of x axis. This transformation will be used below 
in the study of stability of dark solitons. 

Introduction of the parameter c permits one to represent equation ( 63 1 as 

2 



P_ 

>2 



ln(l + 7p) 



1 



(1 + 7)7 



P 



1 



ln(l + 7) 



7(1 + 7) 



c 
2" 



Qip), 



(68) 



where ^ = x + cz. The function Q{p) has a double zero at p = 1 which corresponds to the tails of soliton. Another zero 
at p = pm corresponds to the minimal intensity at the center of soliton, which is, therefore, related to the parameter 
c as 

1/2 



1 



1 



Pn 



2prr. 

7 



li„- 

7 1 



1 + 7 

+ lPn 



1 



1 + 7 



Taking into account Eq. ( 67 ) we find expression for the slope a as a function of p„ 

1/2 



UHl-Pm)'l 



2p„ 



7 l+7Pn7 



1+7 



1 



The slope of the most shallow solitons with p„i — )■ 1 is equal to 



1 



(69) 



(70) 



(71) 



that is it coincides with the Mach cone. 

The profile of the light intensity across the oblique soliton can be obtained by a straightforward numerical integration 
of Eq. ( 68 ) . In Fig. 8 we compare such a profile with the profile of the diffraction pattern obtained by direct numerical 
simulation using original Eqs. (5j6l. Good agreement between these two profiles confirms that the pattern in Fig. 2 
inside the Mach cone indeed consists of oblique dark solitons generated by nonlinear diffraction of the light beam on 
the obstacle. 
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FIG. 8: Profiles of the intensity distributions for x — 100 (daslied line), x = 400 (solid line) and y > obtained from numerical 
solution of the equation ([5| with the nonlinear term given by ([gJi. These profiles are compared with the soliton profiles obtained 
by solutions of Eq. ( |63| l with slope a — 10.58 shown as functions of y at the same values of x (a; = 100 corresponds to "crosses" 
and X = 400 to "circles"). 



V. STABILITY OF OBLIQUE SOLITONS 



The solitons profiles investigated in the preceding Section are reached asymptotically as z — > oo. However, the 
pattern calculated for finite z and shown in Fig. 2 indicate that some oscillations of intensity take place along the 
oblique solitons. Amplitude of these oscillations increases with distance from the obstacle what leads to generation of 
vortices at the end points of solitons. In fact, instability of dark 2D solitons with respect to transverse perturbations 
is well known as well as development of this instability to formation of vortices (see [1]). But in the case of formation 
of dark solitons in the flow of Bose-Einstein condensate past an obstacle it was found [50] that the amplitude of 
oscillations decreases with growth of time at fixed distance from the obstacle for large enough value of the oncoming 
fiow velocity. This suggests that absolute instability of dark solitons transforms into their convective instability in 
the reference frame attached to the obstacle at some critical value of the fiow velocity [31 . This means that wave 
packets built of unstable modes of soliton 's disturbance are convected so fast by the fiow that they cannot develop at 
finite distance from the obstacle. The criterion of transition to the convective instability for Bose-Einstein condensate 
evolving according to the Gross-Pitaevskii equation (see ^) was derived in 31J and here we shall extend the analysis 
of [5P to the photorefractive equation 



A. Shallow solitons (Kadomtsev-Petviashvili approximation) 

The theory is especially simple in the limit of small-amplitude solitons when the GNLS equation (|5| can be reduced 
to the Kadomtsev-Petviashvili (KP) equation by means of standard reductive perturbation theory, which yields 

[-2c,p; + 2c2p^ + (3/'(po) + Pof'iPo)) p'p'i - Ip'ssi] , + clp'yy = 0, (72) 

where p' <^ po denotes the intensity perturbation small compared with the background intensity pQ, x is a, coordinate 
along a soliton and y is a transverse coordinate. We transform it to the standard form by introducing the new variables 



to obtain 



, . x + csz, fi^^, p=-|(3/'(po) + Po/"(po))p' (73) 

ZCg Cc 



(^p, + 3pp^ + = Pfifi- (74) 



As is well known, the KP equation (74 1 has the soliton solution 

s 



cosh^[Vs(^ - sS)] cosli [y^(5 + {cs - s/{2cs))z)] 



(75) 
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where the parameter s is small, 



« 1, 



(76) 



in accordance with the condition that the soliton is shallow. This solution solution is written in the reference frame 
with u — > as x ^ cx). 

The soliton solution (75) is unstable with respect to transverse perturbations |37l |38]. If we perturb the solution 

(77) 



( 75 1 along y axis, 



P = PsiO + ^P^ Sp=W{C)exp{Tz + ipy) , ( = i + cy, 
then in linear approximation we obtain equation for W: 



(78) 



This eigenvalue problem was studied in |38[ 139] where the following spectrum for the instability growth rate was 
obtained 



T{p)^{p/V3)\/s-2p/V3. 



(79) 



Thus, in the reference frame with u ^ at a; — > oo the soliton is absolutely unstable. 

However, we are interested in the behavior of the soliton transformed to the reference frame "attached" to the 
obstacle by substitution ( 64 1 : 



s cosh 



1 + a^ 



X — ay + 



s 

2^ 



\/l+~a^~u] z 



(80) 



The relationship between the soliton parameter s and the slope a follows from the condition that the oblique soliton 
solution does not depend on z: 



1 - 



l + a2 



(81) 



where we took into account ( 13 ) and ( 76 ) . After the transformation to the "obstacle" frame we easily get the dispersion 
relation 



Lli = uj{p) = iip + i 



■ P 



2p Ma 

s = , u = U sm t) = , 

V3' ^ il + j)VT+^ 



(82) 



for waves propagating along oblique soliton with the wave number p. The stability of the soliton is determined by the 
asymptotic behavior of the wave packets built from harmonic waves. Due to the term fip in the dispersion relation, 
the wave packets are convected by the flow along the soliton. If they are convected fast enough, then amplitude 
of the unstable disturbance cannot increase at fixed distance from the obstacle and, as a result, the soliton is just 
convectively unstable 31j. As was shown in [31] . for shallow KP solitons the criterion of transition from absolute to 
convective instability reads 



> s. 



Then substitution of Eqs. (81 1 for s and (82 1 for fj. gives at once 

M > 1. 



(83) 
(84) 



Thus, the shallow solitons are convectively unstable for "supersonic" values of transverse wave vector U. 



B. Deep solitons 



Now we consider stability of soliton solutions of the photorefractive equation ^ with dropped external potential: 
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Stability of solitons for the case of 2D NLS equation (7 = 0) was studied in [HP . We shall write the soliton solution 
of Eq. ( 85 1 in the form 



iz 



1 + 7 



where PsiC) is given by the solution of Eq. (68 1 and 



1 



dC \PsiC) 
The disturbed function ■0 can be taken in the form 

Vs(C) + (V'' + #")exp(z0.(C) 



(86) 



(87) 



1 + 7 



Here ip' and V'" depend on y and z as exp(ipy + Tz) Substitution of Eq. (88 1 into (85 1 and linearization with respect 
to ^' and Ip" yields the linear spectral problem 



L2 A 



r 



r 



where 



Ps \di 2ps 



(89) 



(90) 



1 92 



2d^^ ' 2'" " ' ■ 1 + 7 



1 c2 
2p! 



3ps + 
(1 + 7/5.)^ 



(91) 



19^ 1. 2 
2de ^ 2^"^ 



1 + 7 



lc2 



l + 7Ps 



(92) 



The function ps is considered here as known for a given value of the soliton velocity c; hence the system ( 89 1 can be 
solved numerically which yields the spectrum of the growth rate T — T{p) for all values of c. 

Again, we transform this solution to the reference frame attached to the obstacle and arrive at the dispersion 
relation 



uj{p) = p.p+iT{p). 



(93) 



This equation determines implicitly the function p = p{lo). The type of stability is determined by the location of 
branching points p^r of this function (see, e.g., [H]) where du/dp — what gives the equation 



dV 
dp 



(94) 



which determines the branching point phr as a function of yit at a given value of c. As was shown in |31j . the critical 
value Per of transition from absolute instability to convective one is determined by the condition that the function 
ptr{p) has a branching point at /i = pcr- This gives the equation 



dp'^ 



= 0, 



(95) 



P=Pa 



solution of which gives the critical value Pcr for a given c. Example of the plot of the absolute value of the function 
T(p) is shown in Fig. 9 for 7 = 0.1 and soliton velocity c corresponding to the minimal intensity p,„ = 0.2 and 
calculated by means of Eq. (69). It has an inflection point at p = Pcr in the region where r(p) is purely imaginary; 



thus. Per can be calculated for a set of values of /9,„ in the interval < pm < 1- 
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FIG. 9: Absolute value of the growth rate F as a function of the wave number p of a harmonic transverse perturbation for 
7 = 1 and pm = 0.2. 



When Per is found as a function of pm, we can substitute its value into Eq. (94) to obtain the critical value of p, 
again, as a function of Pm- 



dp 



P=Pa 



Now we substitute the relation 



M 



(i + 7)\/rT^ 

into /i = Ma/ ((1 + 7) vT-H?) to find p, ~ ca which gives the slope parameter as a function of p„ 

Per {Pm) 



a{pm) 



c{pm) 



(96) 



(97) 



(98) 



where c{pm) is given by Eq. (69). At last, substitution of this function into Eq. (97) yields M^r as a function of p,„: 



Mcr{p,n) = (1 +7)c(Pm)Vl +a2(p„). 



(99) 



Equations ( 98 ) and ( 99 ) determine the critical value of Mach number as a function of the slope a in a parametric 



form with < pm ^ 1 playing a role of the parameter. Results of numerical computation of this function for several 
values of 7 are shown in Fig. 10. Below these curves oblique solitons are absolutely unstable and cannot be created 
by the "flow" of light past an obstacle: perturbation of the "flow" behind the obstacle decays into vortices without 
formation of solitons. Above these curves oblique solitons become just convectively unstable and their length grows 
up faster than they decay into vortices. Hence, vortices exist at the end points of solitons only and there is a region 
where the soliton proflle is close to the stationary solution found in the preceding Section which was confirmed by 
numerical simulations. 



VI. CONCLUSION 



In this paper, we have developed the theory of formation of the wave pattern of light propagating through nonlinear 
photorefractive medium with a reflecting wire embedded in the medium. The light beam is supposed to be tilted 
with respect to the wire which creates the "flow of light past an obstacle" analogous to that realized in experiments 
on superfluid flow of Bose-Einstein condensate past an obstacle. An analogy between propagation of light beams and 
superfluid dynamics suggests that diffraction pattern similar to what was observed and predicted theoretically can 
be found in optical experiments. We have shown that the diffraction pattern consists of two regions separated by 
the "Mach cone" outside which the so-called "ship waves" are located while outside this "Mach cone" the nonlinear 
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FIG. 10: Boundary between regions of absolute and convective instabilities for several values of the saturation parameter 7. 



dispersive shocks generating obhque sohton trains are situated. The simplest case when just a single soliton is 
generated is studied in detail. The main parameters of the oblique optical soliton are determined and it is shown that 
it is actually stable (more precisely, convectively unstable) with respect to small transverse perturbations for large 
enough values of transverse wave vector of the light beam. Detailed theory of "ship waves" is also given. All our 
findings are confirmed by numerical simulations. Since optical experiments seem more feasible than the experiments 
with ultra cold gases, one may hope that our predictions could be verified experimentally. 
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